%   Test the following QSEPACK Template Library routines:
%
%            initialize_metric_QSE
%              compute_product_QSE
%      compute_explicit_metric_QSE
%
%   Both validity and speed are tested.  The test takes on the order
%   of one minute.  Printed comments explain the tests being run.

% ******************************************************************************
% Copyright (C) 2008 John Arthur Sidles.  This copyright is assigned  
% to the Institute for Soldier Healing of Seattle, Washington, USA.
% ------------------------------------------------------------------------------
% This source file is part of the QSEPACK Template Library (release 1.0d).  
% The QSEPACK Template Library is free software: you can redistribute it 
% and/or modify it under the terms of the GNU General Public License 
% as published by the Free Software Foundation, either version 3 of 
% the License, or (at your option) any later version. 
% 
% The QSEPACK Template Library is distributed in the hope that it will  
% be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
% GNU General Public License for more details. 
% 
% You should have received a copy of the GNU General Public License along  
% with the QSE Template Library. If not, see <http://www.gnu.org/licenses/>.
% 
% The author, John A. Sidles, can be contacted through the Institute for 
% Soldier Healing, 4532 48th Av NE, Seattle, Washington, 98105 USA, or
% by email at "jsidles@soldierhealing.org"
% ------------------------------------------------------------------------------
% The QSEPACK template library is maintained using the open-source
% literate programming tool nuweb (v1.1.1), by Preston Briggs.
% ------------------------------------------------------------------------------
% The QSEPACK template library is distributed as a gzipped tarball file 
% "QSEPACK_v1_0d.tar.gz" that includes the master documentation file 
% "QSEPACK_v1_0d.pdf". 
% ******************************************************************************

% ------------------------------------------------------------------------------
% scrap: validity test of test_compute_product_QSE.m

clear all; % clear everything (including variables, breakpoints, functions).
clear clc; % clear the screen too.

disp(sprintf( ...
  '\n--------------------------------------------------------------------------------'));
disp(sprintf('Testing the following QSEPACK Template Library routines:\n'));
disp(sprintf('         initialize_metric_QSE'));
disp(sprintf('           compute_product_QSE'));
disp(sprintf('   compute_explicit_metric_QSE'));

rankIn = 5;                    % the rank of the test matrix
dimsIn = [2,3,5,4,3,2,3,2]';   % the dimensionality of the spin states


% ------------------------------------------------------------------------------
% scrap: compute a random input xi-state
xiTest = cell(rankIn,length(dimsIn));
for in=1:rankIn
  for jn=1:length(dimsIn)
    xiTest{in,jn} = complex( ...
      randn([1,dimsIn(jn)]),randn([1,dimsIn(jn)]) ...
      )./ sqrt(2*dimsIn(jn));  % (normalized)
  end
end

% ------------------------------------------------------------------------------
% scrap: compute the structured representation of G
Gtest = initialize_metric_QSE(xiTest);

% ------------------------------------------------------------------------------
% scrap: compute a random test vector
vInTest = complex( ...
      randn([Gtest.rank*Gtest.ncoef,1]), ...
      randn([Gtest.rank*Gtest.ncoef,1]) ...
      );

% ------------------------------------------------------------------------------
% scrap: perform the matrix-vector multiplication
vOutTest = compute_product_QSE(Gtest,vInTest);

% ------------------------------------------------------------------------------
% scrap: evaluate whether the result is valid 
errMag = max(abs(compute_explicit_metric_QSE(Gtest)*vInTest - vOutTest));
disp(sprintf( ...
'\n--------------------------------------------------------------------------------'));
if errMag < 1.0e-10
  disp(sprintf( ...
  'VALIDATED that the following QSEPACK Template Library routines are working:\n'));
  disp(sprintf('         initialize_metric_QSE'));
  disp(sprintf('           compute_product_QSE'));
  disp(sprintf('   compute_explicit_metric_QSE'));
else
  error('QSEPACK:test_compute_product_QSE:validationFailed', '%s', ...
    'The validation of test_compute_product_QSE failed.');
end

% ------------------------------------------------------------------------------
% scrap: speed test of test_compute_product_QSE.m
disp(sprintf( ...
  '\n--------------------------------------------------------------------------------'));
disp(sprintf( ...
  'Testing the speed of the following QSEPACK Template Library routines:\n'));
disp(sprintf('         initialize_metric_QSE'));
disp(sprintf('           compute_product_QSE\n'));

disp(sprintf( ...
  'Preparing the test case described in "Five Promises from Geometric'));
disp(sprintf( ...
  'Quantum Mechanics for Efficient Quantum Simulation" ...\n'));
disp(sprintf( ...
  'Note: the following four timing tests require on the order of one gigabyte'));
disp(sprintf( ...
  '      of free memory, and take on the order of one minute to execute.'));
disp(sprintf( ...
  '\n      ... so be patient ...'));

% ------------------------------------------------------------------------------
% to change the size of the test case, adjust rankIn and dimsIn 
rankIn = 100;               % the rank of the test matrix
dimsIn = 2*ones([500,1]);   % the dimensionality of the spin states
% ------------------------------------------------------------------------------

disp(sprintf( ...
  '\nTiming for a %d X %d (structured) matrix-vector product:', ...
  rankIn*sum(dimsIn),rankIn*sum(dimsIn)));
ticTime = cputime;

% ------------------------------------------------------------------------------
% scrap: compute a random input xi-state
xiTest = cell(rankIn,length(dimsIn));
for in=1:rankIn
  for jn=1:length(dimsIn)
    xiTest{in,jn} = complex( ...
      randn([1,dimsIn(jn)]),randn([1,dimsIn(jn)]) ...
      )./ sqrt(2*dimsIn(jn));  % (normalized)
  end
end
disp(sprintf( ...
  '%8.2f sec to compute a random input xi-state',cputime-ticTime));
ticTime = cputime;

% ------------------------------------------------------------------------------
% scrap: compute the structured representation of G
Gtest = initialize_metric_QSE(xiTest);
timeToComputeG = cputime-ticTime;
disp(sprintf( ...
  '%8.2f sec to compute the structured representation of G',timeToComputeG));
ticTime = cputime; 

% ------------------------------------------------------------------------------
% scrap: compute a random test vector
vInTest = complex( ...
      randn([Gtest.rank*Gtest.ncoef,1]), ...
      randn([Gtest.rank*Gtest.ncoef,1]) ...
      );
disp(sprintf( ...
  '%8.2f sec to compute a random test vector',cputime-ticTime));
ticTime = cputime;

% ------------------------------------------------------------------------------
% scrap: perform the matrix-vector multiplication
vOutTest = compute_product_QSE(Gtest,vInTest);
timeToComputeProduct = cputime-ticTime;
disp(sprintf( ...
  '%8.2f sec to perform a %d X %d (complex) matrix-vector product', ...
  timeToComputeProduct,rankIn*sum(dimsIn),rankIn*sum(dimsIn)));

disp(sprintf( ...
  '\nEstimating time to compute matrix-vector product directly'));
disp(sprintf( ...
  '(we have to estimate, rather than measure directly, because '));
disp(sprintf( ...
  'the explicit %d X %d matrix is too large to hold in memory).', ...
  rankIn*sum(dimsIn),rankIn*sum(dimsIn)));

nDirect = min([rankIn*sum(dimsIn),2500]); % the exact ceiling is not critical
matrixDirect = complex(randn([nDirect,1]),randn([nDirect,1]))* ...
  complex(randn([1,nDirect]),randn([1,nDirect]));
vectorInDirect  = complex(randn([nDirect,1]),randn([nDirect,1]));
vectorOutDirect = complex(zeros([nDirect,1]),zeros([nDirect,1]));
ticTime = cputime;
vectorOutDirect = matrixDirect * vectorInDirect;
timeToComputeReferenceProductDirect = cputime-ticTime;
disp(sprintf(...
  '%8.2f sec for a direct %d X %d (complex) matrix-vector product', ...
  timeToComputeReferenceProductDirect,numel(vectorInDirect),numel(vectorInDirect)));
timeToComputeProductDirectEstimated = timeToComputeReferenceProductDirect * ...
  (rankIn*sum(dimsIn)/numel(vectorInDirect)).^2;

disp(sprintf( ...
  '\nExtrapolating above timing to full-size matrices yields:'));
disp(sprintf( ...
  '%8.2f sec (extrapolated) for a direct %d X %d (complex) product', ...
  timeToComputeProductDirectEstimated,rankIn*sum(dimsIn),rankIn*sum(dimsIn)));

speedRatio = timeToComputeProductDirectEstimated/timeToComputeProduct;
disp(sprintf( ...
  '\nThe gain in speed in computing the matrix-vector product is a factor of %.1f', ...
  speedRatio));

% ------------------------------------------------------------------------------
% scrap: storage test of test_compute_product_QSE.m
disp(sprintf( ...
  '\n--------------------------------------------------------------------------------'));
disp(sprintf( ...
  'Testing the storage of the following QSEPACK Template Library routines:\n'));
disp(sprintf('         initialize_metric_QSE'));

s = whos('Gtest'); numBytesInG = s.bytes; clear s ;
disp(sprintf( ...
  '\nThe number of bytes used to store the %d X %d metric G is %.3g', ...
  rankIn*sum(dimsIn),rankIn*sum(dimsIn),numBytesInG));

testC = complex(1.0,1.0); s = whos('testC'); numBytesInComplex = s.bytes; clear s ;
numBytesInGDirect = numBytesInComplex*rankIn*sum(dimsIn)*rankIn*sum(dimsIn);
disp(sprintf( ...
  'The number of bytes required for direct storage is %.3g', ...
  numBytesInGDirect));
  
storageRatio = numBytesInGDirect/numBytesInG;
disp(sprintf( ...
  '\nThe gain in storage efficiency of the state-space metric is a factor of %.1f', ...
  storageRatio));

% ------------------------------------------------------------------------------
% scrap: validation summary of the QSEPACK Template Library
disp(sprintf( ...
  '\n********************************************************************************'));
disp(sprintf( ...
  'Summary of validation and testing of the QSEPACK 1.0d Template Library '));
disp(sprintf( ...
  'for the environment\"MATLAB %s\" on a \"%s\"',version,computer));
disp(sprintf( ...
  '  * The QSEPACK 1.0d routines appear to be operating correctly.', ...
  speedRatio));
disp(sprintf( ...
  '  * The test case is a matrix-vector product for a %d X %d metric', ...
  rankIn*sum(dimsIn),rankIn*sum(dimsIn)));
disp(sprintf( ...
  '    - The gain in matrix-vector product speed is a factor of %.1f.', ...
  speedRatio));
disp(sprintf( ...
  '    - The gain in metric storage efficiency is a factor of %.1f.', ...
  storageRatio));
disp(sprintf( ...
  '********************************************************************************'));